Nonlocal interactions prevent collapse in Bose-Einstein gases 
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We study the effect of nonlocality on the collapse properties of a self-focusing Nonlinear Schrodinger 
system related to Bose-Einstein condensation problems. Using a combination of moment techniques, 
time dependent variational methods and numerical simulations we present evidences in support of 
the hypothesis that nonlocal attractively interacting condensates cannot collapse when the dominant 
interaction term is due to finite range interactions. Instead there apppear oscillations of the wave 
packet with a localized component whose size is of the order of the range of interactions. We discuss 
the implications of the results to collapse phenomena in negative scattering length Bose-Einstein 
condensates. 
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I. INTRODUCTION 

The problem of collapses in nonlinear wave equa- 
tions has been studied extensively from the mathematical 
point of view ffll. In the collapse phenomenon the am- 
plitude of a physical quantity becomes infinite at a par- 
ticular time and usually the mathematical model is no 
more representative of the physics of the problem. This 
is the case in many particular instances of collapses in 
plasma physics 0], nonlinear optics ||, and many other 
examples. 

In the context of Bose-Einstein condensation there has 
been recent interest on the analysis of collapse prob- 
lems. Condensates near T — OK are modeled by a multi- 
dimensional Nonlinear Schrodinger equation called the 
Gross-Pitaevskii (GP) equation [Q. In this framework, 
the interaction between the constitutive bosons inside the 
condensate is defined in terms of the ground state scat- 
tering length a. When a > the interaction between 
the particles in the condensate is repulsive which reflects 
into a positive self-interaction term in the GP equation, 
whereas for a < 0, the interaction is attractive and self- 
interaction leads to collapse in dimension higher or equal 
than two |^-^]. It is thought that the practical impli- 
cation is that, in this case, the realization of the con- 
densate is limited by a critical number of particles ||] 
above which the condensate is unstable and destroyed by 
the collapse phenomenon 0] . This is the reason why 
most BEC experiments use gases with positive scattering 
length since negative scattering length condensates seem 
to be limited to a very small number of particles. 

In spite of collapse, negative scattering length conden- 



sates have some peculiarities which make them interest- 
ing. One of them is that if the trap is removed in one 
direction, the attractive interaction yields (for a particu- 
lar number of particles) to a self-confined stationary state 
as it has been shown in a previous work Q . In that case, 
the cloud as a whole behaves as a particle and can be 
controlled by acting on it with external fields. These 
solutions can be properly called solitons since they ap- 
pear as a consequence of self-interaction and not merely 
because of the effect of an external potential. These soli- 
tonic condensates would thus be more controllable than 
usual repulsive condensates. 

There is another reason why one would be interested 
in those condensates. Although it is thought that these 
condensates are unstable, this fact has not been really 
tested experimentally. What Rice group experiments es- 
tablished H was that it is not possible to generate a large 
enough negative scattering length condensate. This prob- 
lem was analyzed in more depth in Ref. |io[ were the 
authors used the two fluids model to analyze the con- 
densation process. They found a natural bound on the 
number of particles attainable by the evaporative cooling 
method which is used to grow the condensate. 

However, condensate growth is not the only way to gen- 
erate large negative scattering length condensates. Re- 
cent experimental results jllj show that it is possible to 
use Feshbach resonances to continuously detune the value 
of a from positive to negative values, by means of exter- 
nal magnetic fields. This provides new interest to the 
analysis of attractively interacting condensates since it 
would provide an alternative way to generate them. An- 
other way these large condensates could be obtained is 
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by overlapping two subcritical negative scattering length 
stable Rb or Cs condensates as proposed in Ref. fl2fl . In 
that paper it has been shown that fringe systems gener- 
ated by interference can be used to stabilize large over- 
lapping negative scattering length condensates. However, 
no matter how one generates the supercritical negative 
scattering length condensate the question is: what would 
be its subsequent dynamics? 

To make a complete analysis of that point, i.e. dynam- 
ics of an already generated supercritical condensate near 
T = OK; one should take into account at least the effect 
of new terms correcting the usual Gross-Pitaevskii equa- 
tion. Although there is the possibility of considering the 
problem in more general frameworks such as the two flu- 
ids model or even the full quantum problem we will con- 
centrate here on simpler formalisms which are already 
very rich as will be shown later. Between the possible 
corrections to GP equation are those due to four-particle 
collisions, which were considered in Rcfs. |l3|,|l4|]. How- 
ever the real relevance of those terms is doubtful since 
there is no idea on what could be the value of the coeffi- 
cient of the higher order correction to the cubic nonlinear 
term. Moreover it can be seen that they must be very 
small |Q and probably losses terms are quite relevant be- 
fore this term becomes appreciable. Finally it is not clear 
that one could perform a Taylor expansion of the non- 
linear interaction functional around the minimum point 
in that situation, a requisite which must be satisfied if 
the true interaction potential is to be approximated by a 
polynomial on 

On the other hand, there is a characteristic of the con- 
densate self-interaction which is neglected in the usual 
treatment: the fact that the true interaction potentials 
are nonlocal. This is the way the self-interaction appears 
in the mean field formulation of the Hartree-Fock prob- 
lem for the Boson system [QQ]. In this paper we study 
the effect of nonlocal interaction terms on the collapsing 
properties of the Gross-Pitaevskii equation. We do not 
claim that this will be the most relavant physical contri- 
bution in all the possible situations but we believe that 
it could be in some of them due to the many different 
atomic species which are known to possess negative scat- 
tering length. It is thus interesting to understand what 
is the isolated effect of nonlocal nonlinear terms. Our 
interest will be twofold: firstly from the point of view 
of nonlinear science we are interested in understanding 
what is the effect of such term on the collapse problem; 
secondly, we will consider the problem from the viewpoint 
of Bose-Einstein condensation phenomena. 

The study of nonlocal nonlinear wave equations is non- 
trivial and only recently there have been several studies 
considering particular cases in the contexts of nonlinear 
optics |l6|-[n|, Quantum Mechanics |^0| and other non- 
local wave equations Presently only few rigorous 
results and some qualitative estimates are available for 
very particular cases. To our knowledge, the only related 
field where these nonlocal Nonlinear equations (in par- 
ticular of Nonlinear Schrodinger type) have been studied 



in detail and with mathematical rigor is scattering the- 
ory, but the results are completely disjoint with our pur- 
pose here p^ |. Finally, in the context of Bose-Einstein 
condensation there has been only a work known to us, 
in which non locality effects are considered [p3| . The 
authors use the gradient expansion to obtain a local in- 
teraction energy and then perform a classical time in- 
dependent variational analysis with Gaussian ansatz to 
estimate the possibility that nonlocality enhances stabil- 
ity conditions. An important conclusion made by the au- 
thors is that nonlocal interactions could prevent collapse. 
However, the results obtained so far are based on qualita- 
tive arguments (the gradient expansion is not applicable 
in the region where the narrowing of the wavepacket is 
stopped) and leave open the question about the dynamics 
of the nonlocal interacting condensate. In our paper we 
complement static results with mathematically rigorous 
analysis of the wavepacket width evolution, time depen- 
dent variational analysis, numerical simulations and the 
analysis of the strongly nonlocal limit. Using analysis of 
the dynamics on time we can say on more firm ground 
that nonlocal interactions of a rather general (but prob- 
ably not any) type prevent collapse. 

Our detailed plan is as follows. In Sec. [n| we present 
the model equations with the nonlocal term. A formal 
analysis of the strongly nonlocal limit analysis and wave 
packet width evolution using the moment method is done 
in Sec. III. In Sec. IV we present some analytical ap- 



proximations to obtain a qualitative description of the 
collapse dynamics in all the regimes. In Sec. [V] we present 
the results of numerical simulations of the full model and 
compare the results with analytical predictions of pre- 
vious sections. Finally, in Sec. 
conclusions. 
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we summarize our 



II. MATHEMATICAL MODEL 

The usual theoretical model used to describe a system 
of 2-body interacting bosons of mass m, with a fixed 
mean number TV, trapped in a parabolic potential V(r) 
is the Gross-Pitaevskii (GP) equation Q] 

dt 2m w 



U 



J K{r-r')\y(r')\ 2 dr' 



(1) 



where U = Anhfa/m characterizes the 2-body interac- 
tion and has been extracted from the kernel K for later 
convenience. The normalization for \& is N = J | \T/ 1 2 d 3 r, 
and the trapping parabolic potential is given by V(r) = 
\mv 2 (X 2 x 2 + X 2 y y 2 + X 2 z 2 ) , as discussed e.g. in Ref. [||. 
Here we have kept the nonlocal interaction as it appears 
in the Hartree-Fock theory. Usually the kernel function is 
taken to be a delta function because of the low energy of 
the interaction process. One can easily recover that case 



2 



by imposing that K(r — r') = S(r — r'). In our treatment 
we will keep the nonlocal nature of the kernel and show 
that it gives rise to many new effects. 

Let us write the equations in a natural set of units 
which for our problem is built up from two scales: the 
trap size (measured by the width of the linear ground 
state), ao = yjTi/mv, and its period, \jv. The new 
variables are (xo, yo, zo) = {x,y,z)/aQ,T = ist, V>( r o) = 
\f(r) yJa^/N . With these definitions the equation simpli- 
fies to 



.dih 



u 



J ^(r -r')|V(r')| 2 *' 



(2) 



being U — AirNa/ao. The particular shape of the ker- 
nel function could depend strongly on the energy of the 
interaction and the geometry of the molecules involved 
and probably it is very difficult to know its precise shape 
except for the case of atomic hydrogen, in which BEC 
has been recently achieved |p4j| . 

In this paper we are interested in basic general qual- 
itative results and we do not try to model the specific 
details of the interaction for any particular atom. This 
is why we will concentrate in the simplest case where the 
kernel depends on one parameter e, related to the ker- 
nel "size" which charaacterizes the range of interactions, 
such that 



limine (r) = SM 



(3) 



Moreover, we will analyze molecules (and thus the inter- 
action) which are spherically symmetric. The main im- 
plication of this fact is that the kernel depends only on 
the distance between two atoms, K(\r — r'|). The ques- 
tion we will address in what follows is the possibility of 
blow-up with regular kernels K e . This is a very intricate 
mathematical problem, which even in the simplest 5-case 
(cubic nonlinearity) is not completely solved, and only 
some estimates on collapse conditions exist. Thus, we 
will not try to provide complete rigorous proofs but only 
to join numerical simulations and approximate analytical 
techniques to understand the problem. 



III. STRONGLY NONLOCAL LIMIT 

Let us start with some physical arguments showing 
that the nonlocality can prevent collapse. Although it 
is not essential the algebra is simplified by assuming 
cylindrically (A x = X y = 1) or spherically symmetric 
(A x = A y = \ z = 1) traps. We notice that there exist at 
least two integrals of motion of (0): 

N = I \tb\ 2 dr, (4a) 
11 = I {\Vib\ 2 + J(r)\rb\ 2 + p 2 \^\ 2 ) dr, (4b) 



where p = |r| and 

J(r) = U J K{r-v')\^{v')\ 2 dv'. 



(5) 



N is the number of particles because of the normalization 
imposed on ip and H is associated with the energy of the 
condensate. 

Let us now consider the situation "near" the collapse 
assuming that the wave function is strongly localized. 
More precisely, it will be assumed that the localization 
region of the wave function £ is much less than the range 
of the nonlocal interactions e, £ <C e |2f|. Also it wwill 
be assumed that the kernel funnction is nonsingular at 
r — > 0. Then in the leading approximation one can ap- 
proximate 

J(r) = UNK(r) + 0(£/e) 

Making the natural supposition that if the collapse occurs 
it happens at r = (i.e. at the minimum of the confining 
potential) one finds the following linear equation for the 
wave function 



.dih 



1 



V 2 Q ib + UNK{Q)ib 



(6) 



In Eq. (|g) we have neglected the confining potential, 
since it is of order of l 2 v <§; 1. 

As it is evident, Eq. (||) does not display collapse. 
Moreover the dispersion will lead to spreading out of any 
initially localized wave packet. 

On the other hand if the condensate width is much big- 
ger than the range of interactions £ 3> e, the ^-function 
limit is applicable and the dominant nonlinearity will lead 
to collapse. Thus we conclude that there must exist two 
different tendences: spreading for I -C e and narrow- 
ing for f > £ of the wave packet. This must lead to 
oscillations of the wave packet width between different 
scales, the smaller one being of the order of the inter- 
action range. As a matter of fact the limit f < e is 
just opposite to the limit of local interactions and can 
be interpreted as the case of infinite range interactions, 
K(r) = K(0). So we see that there is a scale at which 
the interaction is attractive and a local scale at which it 
is repulsive. It seems natural that this will be reflected 
on oscillation in the wavefunction dynamics, a fact which 
we will analyze later. 

It is also possible to argue that the collapse does not 
occur in this system using the more conventional lan- 
guage of momenta (0,^6|- To this end let us define the 
mean squared width of the wave packet 



(p 2 )=2vr(n-l) / \^(p)\ 2 p n+1 dp 



(7) 



where n = 2, 3 is the spatial dimension. Then it is a 
straightforward algebra to obtain 
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rfr 2 



(p 2 ) = 2H- 4A 2 ( ( o 2 ) - 47r(n - l)U 



J(p) 



1 



pJ'(p) 



p n ' X dp (8) 



In the limit of infinite range of interactions (or very lo- 
calized wavefunctions) one finds 



— (p 2 ) = 2H-4X 2 (p 2 



1 



n(n — 1) 
The solution of Eq. (§) reads 

(p 2 ) = pi sin(4T + 0q) + C 



UN 2 



(9) 



(10) 



where C = 2H + UN 2 



Po and </>o are real con- 



7r(n^Ty' 

stants. It follows from the energy conservation law that 



C > 2(p 2 ) 



(11) 



Hence the minimal value of (p 2 ), given by C— p\ is always 
positive (p 2 ) > since it follows from ( |Tl| ) that C/4 > 
po. This result rules out the possibility that it could 
exist a collapse in which all the wavefunction became 
concentrated on one particular point (e.g. a delta-like 
singularity) . 



IV. ANALYTICAL RESULTS FOR THE 
GENERAL CASE 

A. Exact results for the center of mass 



It has been pointed out in previous works 
the center of mass defined by 



(r > = J i#| 2 dr 



28 that 



(12) 



performs harmonic oscillations no matter what the non- 
linear interaction is. This result is also valid for the non- 
local interaction case under very general conditions as 
follows from Erhenfest theorem extended to this case so 
that 



— r < x > +X < x >= 0, 
dr z 

d 2 , 2 
^-2 < y > +\ 2 < y >= 0, 

d 2 
dr 2 



—z < z > +X 2 Z < z >= 0. 



(13a) 
(13b) 
(13c) 



To obtain more exact information on this problem a pos- 
sibility is to use the moment method in radial 
symmetry. Moreover it can be seen that it does not pro- 
vide exact results for our case. Even the use of the uni- 
form divergence approximation p9| ] is not possible and 
the only possibility is to restrict to the classical time de- 
pendent variational techniques, which we consider in the 
following subsection. 



B. Time dependent variational formalism 

Following the standard procedure we first identify a 
Lagrangian density for problem (Q), which is 



C 



ih / d** 
— \p 

2 V dt 



~dt 



2m 



J A'(r - r / )|^ r (r)| 2 |^(r')| 2 rfr', (14) 



where the asterisk denotes complex conjugation. That 
is, instead of working with the Gross-Pitaevskii equation 
we can treat the action, 



S 



J Cd 3 rdt = jf 



L(t)dt, 



(15) 



where ti and t / are initial and final moments of time, and 
study its invariance properties and extrema, which are in 
turn solutions of Eq. (|l|) . 

To simplify the problem, we restrict the shape of the 
function f to a convenient family of trial functions and 
study the time evolution of the parameters that define 
that family. A natural choice, which corresponds to the 
exact solution in the linear limit (U — 0) and provided 
quite good results in our previous works |^,^| is a n- 
dimensional Gaussian-like function 



-lv - Vcm} 2 
2w 2 



irjar, + ir] 2 (3 v 



(16) 



where A (amplitude), w v (width), a v (slope - speed), 
/3ri (square root of the curvature radius) and 770 (center 
of the cloud) are free parameters. We are considering n 
as a free parameter that can be set to two or three de- 
pending on the dimensionality of the problem considered. 
BEC systems are in general three dimensional; however, 
in certain situations a two dimensional condensate can 
be considered a good theoretical model |5l| and is easier 
to compare with numerical simulations of Eq. (|l]) . 

The procedure for deriving equations for the parame- 
ters has been described in previous works J6p7| and will 
not be repeated here. 

To go on with the analysis a particular shape must be 
chosen for the kernel function. We will consider a simple 
n-dimensional Gaussian kernel of the form 



K(r) 



( 1 



V27re : 



n/2 



-r 2 /2e 2 



(17) 



The main equations obtained when computing the evolu- 
tion equations for these kind of systems are those related 
to the width. To do so it is useful to introduce a set 
of rescaled variables for time, r = id, and the widths, 
w v = aov v , (77 = x, y, z). For a n-dimensional condensate 
the equations are found to be 



4 



d v k 2 1 

-3-5- + \ k Vk = o" 



n 



8 2 ) 



1/2 > 



(18) 



where P = \j2/nNa/ao (strength of the atom-atom in- 
teraction) and 5 = e/ciQ. The remaining parameters 
a n , f3 v obey separate equations but essentially can be 
computed independetly once the Vk(t) are known. 
The system (Fl8) is generated by the Hamiltonian 



H 



IE 



\ 2 2 



1 



(v 2 +S 2 ) 



1/2' 



(19) 



In terms of the system (|18| ) collapse corresponds to the 
behavior when v k — > 0. As a matter of fact the ex- 
plicit expression of the Hamiltonian already shown that 
collapse is prevented by nonlocal interactions (i.e. by 
nonzero S) in the framework of this simple model. Indeed 
one can easily see that H has a lower bound H > —P/8 3 . 



C. 2D case: Equilibrium points and minimum width 

In the present section we will concentrate on two di- 
mensional condensates for the sake of comparison with 
numerical simulations of Eq. (Q) , which will be presented 
in Sec. [v| In this case the equations are 



Pv x 



d 2 V x 2 1 

+ \V X - 

f^y+x 2 v =1 + ?S 

dT 2 y y V 3 , („2 + £2)3/2 ( W 2 + £2)1/2' 



(20a) 



(20b) 



Let us consider the case when v x = v y — v, which cor- 
responds to a cylindrical symmetry around the center of 
the wavefunction. When the solution has the symmetry 
of the external potential (i.e. t]cm = 0) this corresponds 
to the usual cylindrically symmetric case since the wave- 
function amplitude depends only on r. In this case the 
equations are simpler 



d 2 ^ 

dT 2 



1 



Pv 



{v 2 + 6 2 y 

This equation can be obtained from the potential 

P 



V(v) 



1 



1 



2(v 2 +S 2 ) 



(21) 



(22) 



which is plotted in Figs, [l] and || for some particular 
parameter values. It can be seen that even when P is 
negative, corresponding to negative scattering length, the 
potential is repulsive at the origin so that no blowup is 
possible. 
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FIG. 1. Potential V{v) for P = -10,5 = 0.1 (a) Detail of 
the small scale (b) Large scale 
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FIG. 2. Potential V(v) for P = -10,6= 0.01 (a) Detail of 
the small scale (b) Large scale 

The case of interest for us corresponds to small 5 val- 
ues. In that limit the potential has two scales. For 
v = 0(S) the parabolic term can be neglected since the 
last two terms are dominant. For v = O(oo) — 0{po/5) 
the last term can be neglected (at least when P is not 
large) and the potential is dominated by the parabolic 
term. The existence of two length scales is also clear in 
Figs. @ and |. 

The equilibrium points v* are the solutions of the al- 
gebraic equations 



1 



Pv 



5 2 Y 



(23) 



The 5 = case was analyzed in Ref. [pj; now the situa- 
tion is quite different since collapse is not possible. In a 
generic case one can show that there exists only one pos- 
itive root of Eq. ((H). To this end we introduce a new 
variable z, v/S = e z and rewrite Eq. (B3T) in the form 



sinh(2(z + zq)) cosh 2 z - 



-Pe 



2z 







(24) 



where 5 = e z ° . The left hand side of this equation is a 
monotonic function of z and thus there exists only one 



5 



real root of Eq. (g4j) . Moreover the root is finite for 5 =/= 
what means that the collapse (v — 0) is not possible. 
This result is in agreement with the consideration of Sec. 

in. 

Let us concentrate on the small 5 case, to do so let us 
first define a new variable q — (v/S) 2 > 0, so that Eq. 
(p3|) becomes 



q 2 (q + lf 



(25) 



Since we are going to deal with small <5 values we can 
neglect the left hand side |32| and obtain after some alge- 

provided 



bra the only equilibrium point as g* 
|P| > 1, which implies that 



\p\-l ' 



5 



(26) 



and is consistent with our previous assumptions on the 
linear term. The fact that the equilibrium point depends 
linearly on S is interesting and provides a first estimate 
of the order of magnitude of the turning point of the po- 
tential which is then roughly 0(5) for a wide range of 
initial energies. We insist that this is a lower bound for 
the minimum width since it could happen (as it does) 
that there is only part of the solution of Eq. (|l]) which 
tends to collapse and then the contribution of the non- 
collapsing part will make the width finite. In that case 
our estimate would roughly correspond to the width of 
the collapsing peak. 

If |P| w 1 Eq. ( p6| ) cannot be applied since then the 
denominator can be small. In this case one obtains the 
law?;* = (2 1 / 2 <5) 1 / 3 . 

It is also possible to compute the frequency of the small 
oscillations around the potential minimum, which defines 
a time scale that could be of the order of oscillations 
found in the condensate dynamics due to the competition 
between nonlocal dispersion and trapping forces 



n 




(27) 



An important feature of this formula is that the frequency 
grows as — > 0, i. e. as the range of interactions 5, 
goes to zero. We will check later the validity of these 
predictions. 



D. 3D case: Equilibrium points and minimum width 

Considering again the simplifying assumption that 
v x = v y = v z = v the dynamical equations in this case 
are 



dfv_ 
'dP 2 



1 



Pv 



V* („2 + £2) 



5/2 



(28) 



and the potential 



V(v) 



1 



1 

2^2 



p 



3 (v 2 + S 2 ) 



3/2 



(29) 



As before, collapse can be ruled out since the singular- 
ity has been removed. However the equilibrium point 
satisfies a more complicated equation 



1 



Pv 



{v 2 + S 2 ) 



5/2 



Defining again q — v 2 /S 2 the reduced equation is 



r4 1 p 

5 = s + 

1 5{q+. 



.5/2 



(30) 



(31) 



and the S dependence of q is nontrivial. Without making 
any approximations Eq. (pWj can be written as 



S 2 (q+lf (5Y-1) -PV = 



(32) 



whose solutions for each (P, 6) pair provide the right equi- 
libria v*. Now using the z- variables as in the previous 
subsection one can prove that there exists only one pos- 
itive root of Eq. ( |32| ) . 

It is possible to investigate the orders of the different 
terms in ( |3l| ) and the only simplifying assumption is that 
ijCl, since we expect now that collapse is stronger than 
before as is usual in three-dimensional problems. Using 
this assumption we find that S/q 2 = — P and thus 



<5 5/4 

TpV4 



(33) 



which is smaller than the two dimensional equilibrium 
width. This is an indication that even though it seems 
plausible that collapse does not take place the compres- 
sion of the width due to frustrated collapse process is 
stronger than that of the two-dimensional case. 



V. NUMERICAL RESULTS 



The theoretical analysis of Sects. Ill and IV share 



two common conclusion: (i) There should be a limit on 
the minimum width of the wavepacket, a fact that could 
rule out collapse and (ii) there must exist oscillations 
on the wave packet. In order to test these results and 
the other predictions related to the dynamics of the con- 
densate which we have obtained during our approximate 
variational analysis we have integrated numerically Eq. 
(||). In our numerical simulations we start with Gaus- 
sian initial data and then compute the solution using a 
symmetrized second order in time Fourier pseudospec- 
tral method. Typical grid sizes range from 128 x 128 to 
512x512 points. Simulation times in adimensional units 
where of the order of 50-100 (integration step At = 0.01) 



G 



which correspond to physical values of the order of sec- 
onds, which are about the lifetimes of the condensates 
and dynamically allow span a time interval in which the 
essential features are to be captured. 

In our numerical simulations we studied the region of 
U values contained between the linear case (U = 0) and 
a high value of U = 100, which is about ten times above 
the threshold for collapse in the local case || . In practice 
we worked with S 2 G [0.005, 0.1]. 

The main conclusion of our analysis is that up to the 
precision of the computation we can conclude that there 
is no collapse in all the situations analyzed. When S is 
further reduced or the nonlinear coefficient U is increased 
above a certain limit we cannot decide if collapse is pos- 
sible since then the spatial resolution of the scheme is 
not enough to know whether collapse is real or a mere 
numerical artifact. This is a problem of the numerical 
simulation which is quite difficult to avoid since as will 
be discussed in detail later one needs at the same time a 
large integration region to capture the global wavefunc- 
tion shape and a fine grid to resolve the small amplitude 
peaks thus making the problem computationally hard. 
On the other hand, it is not easy to consider the radi- 
ally symmetric geometry as it is done in local collapse 
problems ||H|||. 

Results of a typical simulation are shown in detail in 
Fig. ||. Two things are clear: firstly that the wave packet 
width oscillates with a minimum width of order 0(1), 
secondly the maximum amplitude of the condensate per- 
forms oscillations with a dominant frequency. 



(P 2 ) L5 




FIG. 3. Width (a) and amplitude (b) oscillations in a su- 
percritical condensate with U = —20, S 2 = 0.05. The frus- 
trated collapse events appear as peaks in the condensate am- 
plitude (maximum height). 



One could expect that, according to Eq. (|26|), the min- 
imum width of the wavepacket were of order S. However 
that would be the case if all the wave function were in- 
volved in the frustrated collapse events, which is not the 
case. As it happens in collapse in the local Nonlinear 
Schrodinger equation only part of the wave function takes 
part on the squeezing dynamics, the remaining part be- 
ing the responsible of the finite size width. In our case 
this is very clear in Fig. || where the two contributions 



to the solution are clearly seen, one being of order 0(1) 
and the other corresponding to a smaller scale t c . Also 
it is interesting how the low amplitude extended ("inert 
part") oscillates according to trap frequency while the 
peak dynamics is ruled out by the nonlocal dynamics, a 
phenomenon which is clear in the width oscillations of 
Fig. |^(a), were the two frequencies are present in the dy- 
namics. In fact, the local oscillations were predicted in 
Sec. [n| and appear also as oscillations in the potential 
well in the variational formalism. They correspond to 
"frustrated collapse" events since if the nonlocality were 
not present the concentration dynamics would not stop 
at scale £ c but would continue up to infinity. 




y 




FIG. 4. (a) Trasversal section of the condensate amplitude 
\ip(0, y)\ during a frustrated collapse event for t ~ 0.4. The 
two contributions to the solution are very clear. The spatial 
spectrum (b) has also signatures of the existence of two scales. 
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FIG. 5. Oscillation frequency of the condensate (circles) 
for U = — 20 and different 8 values against the variational 
estimate (solid line). 

It is remarkable that even when the variational method 
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does not take into account the existence of two scales at 
least some of its predictions match reasonably well with 
the simulated dynamics of Eq. (^) . For example Fig. || 
shows the frequency of the amplitude oscillations of the 
condensate as a function of S (circles). The variational 
estimates (solid line) are in qualitative agreement. This 
result cannot be improved in the framework of simple 
variational analysis since it corresponds to the Gaussian 
ansatz and is derived for the oscillations near the bottom 
of the potential well. 

Other result of the variational analysis which is quite 
relevant is the size of the small scale £ c generated during 
the partial compression of the wave packet, which should 
be of order S for a two dimensional condensate. To check 
it we have estimated the width of the collapsing peak 
using the small scale of the spatial representation of ip 
plot and the long scale of spatial spectra for different 
values of 5. Our results are summarized in Fig. ^, were 
it can be seen how the width of the peak depends linearly 
on S, which again supports the variational analysis. 




ble and no asymmetric instabilities grow when one starts 
with symmetric initial data. In fact, in the context of 
local GP equation collapse phenomenon seems to be an 
essentially radially symmetric process j37],[38| . 

The analytical predictions have been tested with a nu- 
merical scheme. Up to the range in which the numerical 
simulations can be trusted we have not observed collapse 
in the system which supports our theoretical analysis, 
even when the last concentrates on the case when all the 
solution collapses (both the strongly nonlocal limit and 
the variational method). 

For BEC the main implication is that if one were able 
to generate large negative scattering length condensates 
as discussed in the introduction, i. e. by overlapping of 
two subcritical ones or by switching the scattering length 
from positive to negative, one could obtain an oscillating 
but non collapsing condensate provided there were no 
additional physical effects to be taken into account. In 
fact since it seems that the Hamiltonian is bounded below 
one could construct a ground state for this system so that 
stationary states can probably be obtained. 

From a fundamental point of view our analysis is the 
first systematic study of the behavior of the effect of non- 
local interacting Bose-Einstein condensates taking into 
account different aspects: variational analysis, moment 
approximations and numerical techniques. We hope this 
study will be valuable in the search for stable attractively 
interacting condensates. 
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FIG. 6. Minimum width of the "collapsing" part £ c (width 
of the small scale peak) as a function of S for U = —20. 
The circles represent the values obtained by the numerical 
simulations of Eq. (j2|) and the solid line is a least squares 
interpolant, which proves the linear dependence on 5. 
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VI. CONCLUSIONS 



In this paper we have presented analytical and numer- 
ical evidences that the effect of nonlocal interactions is to 
prevent collapse when this is the dominant physical inter- 
action process in negative scattering length Bose-Einstein 
condensates. The analytical tools combine the exact 
analysis of the highly nonlocal limit, moment analysis 
and collective coordinate analysis. Even though there is 
no a rigorous proof, the predictions of all methods point 
out that collapse is excluded when nonlocal (and nonsin- 
gular) interactions are taken into account. 

Most of the analytical predictions have been based 
upon the assumption of cylindrical or spherical symme- 
try of the wavefunction although it is easy to generalize 
the results to other symmetries. From the numerical so- 
lutions one sees that radially symmetric solutions are sta- 
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